source("load_packages.R")
source("extract.custom.R")
source("functions.R")
source("formulas.R")
source("variable_labels.R")
library(stargazer)
library(remotes)
#remotes::install_github("ChandlerLutz/starpolishr") #install if you haven't
library(starpolishr)
library(dplyr)
library(countrycode)
library(estimatr)
library(MatchIt)
library(RobustSE)

##################### NOTE ####################################
## Only models using stargazer outputs are used in our paper ##
###############################################################


vars_core_set <- all.vars(as.formula(bf[[4]]))

sy <- 1994
ey <- 2013

load("aiddem_datasets.RData")
dX <- d1
ncountries <- dX %>% filter(smpl == 1) %>% select(cc) %>% unique %>% nrow


vdemdata<-read.csv("V-Dem-CY-Core-v12.csv")
vdemdata<-vdemdata%>%filter(year>=1993)
v2<-vdemdata%>%select(v2elsuffrage,year,country_text_id)
dx2<-left_join(dX,v2,by=c("year"="year","cc"="country_text_id"))

## Democracy 7 years later
lag_dem<-vdemdata%>%select(v2x_polyarchy,year,country_text_id)
colnames(lag_dem)<-c("polyarchy_7ahead","year","cc")
lag_dem$year<-lag_dem$year-7
dx2<-left_join(dX,lag_dem,by=c("year"="year","cc"="cc"))
dx2$polyarchy_7ahead<-dx2$polyarchy_7ahead*100
#dx2$v2psplats_5ahead<-dx2$v2psplats_5ahead*100
#dx2$v2psbars_5ahead<-dx2$v2psbars_5ahead*100

r.ols.1 <- vector("list", 2)
r.ols.1[[1]] <- felm(polyarchy_7ahead ~ l.CMgnh |
                       cnamef + periodf | 0 | cnamef
                     , data = dx2 %>% filter(smpl == 1))
r.ols.1[[2]] <- felm(polyarchy_7ahead ~ l.CMgnh +
                       l.CMgap.log + l.CMenh +
                       l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                       cnamef + periodf | 0 | cnamef
                     , data = dx2 %>% filter(smpl == 1))
r.ols.2<-vector("list",2)
r.ols.2[[1]]<- felm(polyarchy_7ahead ~
                      l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf | 0 | cnamef
                    , data = dx2 %>% filter(smpl == 1))
r.ols.2[[2]]<- felm(polyarchy_7ahead ~ IA_l + l.CMgnh +
                      l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf | 0 | cnamef
                    , data = dx2 %>% filter(smpl == 1))


#screenreg(r.ols.1
#          , custom.coef.map =
#            list(
#              "l.CMgnh" = vln[2],
#              "l.CMgap.log" = vla[2],
#              "IA_l" = vlies[2],
#              "l.CMenh" = vln[3],
#              "l.CMeap.log" = vla[3],
#              "l.pop.log.r" = vlctrls[1],
#              "l.gdpcap.log.r" = vlctrls[2],
#              "l.war25" = vlctrls[3]
#            )
#          , ci.force = TRUE
#          , digits = 2
#)

#screenreg(r.ols.2
#          , custom.coef.map =
#            list(
#              "l.CMgnh" = vln[2],
#              "l.CMgap.log" = vla[2],
#              "IA_l" = vlies[2],
#              "l.CMenh" = vln[3],
#              "l.CMeap.log" = vla[3],
#              "l.pop.log.r" = vlctrls[1],
#              "l.gdpcap.log.r" = vlctrls[2],
#              "l.war25" = vlctrls[3]
#            )
#          , ci.force = TRUE
#          , digits = 2
#)

# latex output 
f1 <- felm(polyarchy_7ahead ~ l.CMgnh |
             cnamef + periodf | 0 | cnamef
           , data = dx2 %>% filter(smpl == 1))
f2 <- felm(polyarchy_7ahead ~ l.CMgnh +
             l.CMgap.log + l.CMenh +
             l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
             cnamef + periodf | 0 | cnamef
           , data = dx2 %>% filter(smpl == 1))

f3<- felm(polyarchy_7ahead ~
            l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
            cnamef + periodf | 0 | cnamef
          , data = dx2 %>% filter(smpl == 1))
f4<- felm(polyarchy_7ahead ~ IA_l + l.CMgnh +
            l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
            cnamef + periodf | 0 | cnamef
          , data = dx2 %>% filter(smpl == 1))
stargazer(f1,f2,f3,f4, covariate.labels =c("No. of dem. donors * Dem. aid (m USD; log) (est.)",
                                           "Number of democracy donors",
                                           "Democracy aid (m USD; log)",
                                           "Number of economic donors",
                                           "Economic aid (m USD; log)",
                                           "Pop. (log) x 10",
                                           "GDP p.c. (log) x 10",
                                           "Civil conflict"),
          omit.stat=c("f", "ser")
)

#plot
r1 <- vector("list", 1)

r1[[1]] <- felm(polyarchy_7ahead ~ l.CMgnh + l.pop.log.r + l.gdpcap.log.r + l.war25
                | cnamef + periodf | 0 | cnamef, data = dx2 %>% filter(smpl == 1))
r1[[2]] <- felm(polyarchy_7ahead ~ l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
                | cnamef + periodf | 0 | cnamef, data = dx2 %>% filter(smpl == 1))
r1[[3]] <- felm(polyarchy_7ahead ~ l.CMgnh + l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
                | cnamef + periodf | 0 | cnamef, data = dx2 %>% filter(smpl == 1))
r1[[4]] <- felm(polyarchy_7ahead ~ l.CMgnh : l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25
                | cnamef + periodf | 0 | cnamef, data = dx2 %>% filter(smpl == 1))
r1[[5]] <- felm(polyarchy_7ahead ~ l.CMgnh * l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                  cnamef + periodf | 0 | cnamef, data = dx2 %>% filter(smpl == 1))

## Estimated coefficient of the number of democracy donors on democracy -----
## conditional on the amount of democracy aid

rp <- lm(polyarchy_7ahead ~ l.CMgnh * l.CMgap.log + l.pop.log.r +
           l.gdpcap.log.r + l.war25 + cnamef + periodf,
         data = dx2 %>% filter(smpl == 1))

# cluster se at country level

interplot(rp, var1 = "l.CMgnh", var2 = "l.CMgap.log",
          sims = 5000, plot = FALSE, hist = TRUE, xmin = 0) -> ip
# get ratio of non-clustered to clustered coefs:
summary(r1[[5]])$coefficients[1, 2] / (1 * summary(rp)$coefficients[2, 2]) -> iprate
# correct se by ratio:
ip$ub <- ip$coef + ((ip$ub - ip$coef) * iprate)
ip$lb <- ip$coef - ((ip$coef - ip$lb) * iprate)

names(ip)[1:2] <- c("fake", "coef1")  # bug in interplot(); must rename

# plot
interplot(ip, hist = TRUE, var2_dt = dX$l.CMgap.log[dX$smpl == 1]
) +
  scale_x_continuous(labels = c("0", "10", "100", "1,000", "")) +
  annotation_logticks(sides = "b") +
  xlab('Democracy aid (m USD)') +
  ylab('Coefficient for the number\n of democracy donors') +
  theme_bw() +
  geom_hline(yintercept = 0, linetype = 2)

## Lagged outcome model
#first create lagged outcome 
dx2<-left_join(dX,v2,by=c("year"="year","cc"="country_text_id"))
dx2.l<-dx2
dx2.l<-dx2.l%>%select(cc,year,v2x.polyarchy.n)
dx2.l$year<-dx2.l$year+1
dx2.l<-dx2.l%>%rename(v2x.polyarchy.n_lag1=v2x.polyarchy.n)
dx2<-left_join(dx2,dx2.l,by=c("cc"="cc","year"="year"))
fit2.1<-lm(v2x.polyarchy.n~ v2x.polyarchy.n_lag1+ l.CMgnh, data = dx2)
fit2.2<-lm(v2x.polyarchy.n~ v2x.polyarchy.n_lag1+IA_l + l.CMgnh +
             l.CMgap.log , data = dx2)
fit2<-lm(v2x.polyarchy.n~ v2x.polyarchy.n_lag1+IA_l + l.CMgnh +
           l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25, data = dx2)
fit2<-glm(v2x.polyarchy.n~ v2x.polyarchy.n_lag1+IA_l + l.CMgnh +
            l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25, data = dx2,family=gaussian)


fit2.r<-lm_robust(v2x.polyarchy.n~ v2x.polyarchy.n_lag1+IA_l + l.CMgnh +
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25, data = dx2, clusters = cnamef, se_type = "CR2")
summary(fit2.r)
summary(fit2)  #classical and robust SEs are similar, so there's no model misspecification
fit2.r

stargazer(fit2.1,fit2.2,fit2, covariate.labels =c("Democracy in the previous year",
                                                  "No. of dem. donors * Dem. aid (m USD; log) (est.)",
                                                  "Number of democracy donors",
                                                  "Democracy aid (m USD; log)",
                                                  "Pop. (log) x 10",
                                                  "GDP p.c. (log) x 10",
                                                  "Civil conflict"),
          omit.stat=c("f", "ser")
)

GIM(fit2,full=FALSE)

## who has more donors
dx2.d<-dx2
dx2.d<-dx2.d%>%select(v2x.polyarchy.n_lag1,l.CMgap.log ,l.CMenh,l.CMeap.log,l.pop.log.r ,l.gdpcap.log.r,l.war25,year,cc)
dx2.d$year<-dx2.d$year+1 #lag explanatory variables
colnames(dx2.d)<-c("v2x.polyarchy.n_lag2", "l.CMgap.log_lag1","l.CMenh_lag1","l.CMeap.log_lag1","l.pop.log.r_lag1","l.gdpcap.log.r_lag1","l.war25_lag1","year","cc")
dx2<-left_join(dx2,dx2.d,by=c("cc"="cc","year"="year"))
fit_donor1<-felm(l.CMgnh~ v2x.polyarchy.n_lag2 |
                   cnamef + periodf | 0 | cnamef
                 , data = dx2 %>% filter(smpl == 1))
fit_donor2<-felm(l.CMgnh~ v2x.polyarchy.n_lag2+l.CMgap.log_lag1 + l.CMenh_lag1+l.CMeap.log_lag1 |
                   cnamef + periodf | 0 | cnamef
                 , data = dx2 %>% filter(smpl == 1))
fit_donor3<- felm(l.CMgnh~ v2x.polyarchy.n_lag2+l.CMgap.log_lag1 + l.CMenh_lag1+l.CMeap.log_lag1+l.pop.log.r_lag1 + l.gdpcap.log.r_lag1 + l.war25_lag1 |
                    cnamef + periodf | 0 | cnamef
                  , data = dx2 %>% filter(smpl == 1))

stargazer(fit_donor1,fit_donor2,fit_donor3, covariate.labels =c("Democracy in the previous year",
                                                                "Democracy aid (m USD; log)",
                                                                "Number of economic donors",
                                                                "Economic aid (m USD; log)",
                                                                "Pop. (log) x 10",
                                                                "GDP p.c. (log) x 10",
                                                                "Civil conflict"),
          omit.stat=c("f", "ser")
)

## Matching 

dX3<-dx2%>%
  select(v2x.polyarchy.n,l.CMgnh,l.CMgap.log ,l.pop.log.r ,l.gdpcap.log.r, l.war25,v2x.polyarchy.n_lag1)%>%
  na.omit()
dX3$l.CMgnh_bin<-ifelse(dX3$l.CMgnh>5,1,0)
maha_match1<-matchit(l.CMgnh_bin~ l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25,
                     data=dX3,distance="mahalanobis",replace=T,method="nearest")
m.data1 <- match.data(maha_match1)

maha_match2<-matchit(l.CMgnh_bin~ v2x.polyarchy.n_lag1+l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25,
                     data=dX3,distance="mahalanobis",replace=T,method="nearest")
m.data2 <- match.data(maha_match2)
fit_maha1<-lm(v2x.polyarchy.n~ l.CMgnh_bin +
                l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25,data=m.data1)
summary(fit_maha1)

fit_maha2<-lm(v2x.polyarchy.n~ l.CMgnh_bin +v2x.polyarchy.n_lag1+
                l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25,data=m.data2)
summary(fit_maha2)

stargazer(fit_maha1,fit_maha2, covariate.labels =c("Many donors","Democracy in the previous year",
                                                   "Democracy aid (m USD; log)",
                                                   "Pop. (log) x 10",
                                                   "GDP p.c. (log) x 10",
                                                   "Civil conflict"),
          omit.stat=c("f", "ser")
)

## Conditionality
c.data<-read_rds("prop_tied.RDS")
c.data<-left_join(dX,c.data,by=c("cc"="cc","period"="period"))
f_c<-felm(v2x.polyarchy.n ~l.CMgnh + l.CMgap.log +l.CMenh+l.CMeap.log+
            l.pop.log.r + l.gdpcap.log.r + l.war25 +prop_tied |
            cnamef + periodf | 0 | cnamef
          , data = c.data %>% filter(smpl == 1))
summary(f_c)
stargazer(f_c, covariate.labels =c("Number of democracy donors",
                                   "Democracy aid (m USD; log)",
                                   "Number of economic donors",
                                   "Economic aid (m USD; log)",
                                   "Pop. (log) x 10",
                                   "GDP p.c. (log) x 10",
                                   "Civil conflict","Prop. Tied Aid"),
          omit.stat=c("f", "ser")
)